Fit density models to cod of different sizes and calculate spatial overlap with prey
Author
Max Lindmark
Published
August 26, 2023
Load packages & source functions
Code
# Using dev branch of sdmTMB# remotes::install_github("pbs-assess/sdmTMB", dependencies = TRUE)# with_libpaths(new = "/Users/maxlindmark/Dropbox/Max work/R/sdmTMB-zeta-intercept/",# install_github("pbs-assess/sdmTMB", ref = "zeta-intercept"))library(sdmTMB, lib.loc ="/Users/maxlindmark/Dropbox/Max work/R/sdmTMB-zeta-intercept")pkgs <-c("tidyverse", "RCurl", "viridis", "devtools", "terra", "readxl", "mapplots", "tidylog") # minpack.lm needed if using nlsLM()if(length(setdiff(pkgs,rownames(installed.packages()))) >0){install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T) }invisible(lapply(pkgs, library, character.only = T))# Packages not on CRAN or dev version# remotes::install_github("pbs-assess/sdmTMB", dependencies = TRUE)library(sdmTMB)# devtools::install_github("seananderson/ggsidekick") # not on CRAN library(ggsidekick)theme_set(theme_sleek(base_size =20))# Source code for map plots# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)# Source code for map plots# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")options(ggplot2.continuous.colour ="viridis")# Source overlap functionsdevtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/overlap-indices.R")# Set pathhome <- here::here()
Read spatiotemporal predictions from 01-fit-sdm.qmd
# Use only one of them because they are so similarcod_pred <-read_csv(paste0(home, "/output/pred_cod.csv"))
Rows: 943428 Columns: 38
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): substrate, month_year, ices_rect
dbl (35): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_d...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# TODO: Update years!# Read data on rectangle levelspr <-read_xlsx(paste0(home, "/data/BIAS/N and B per Rect. 1991-2020.xlsx"),sheet =4) |>mutate(sub_div =ifelse(Sub_Div =="28_2", "28", Sub_Div)) |>filter(sub_div %in%c("24", "25", "26", "27", "28")) |>rename("ices_rect"="RECT","Year"="ANNUS") |>mutate_at(vars(`1`, `2`, `3`, `4`, `5`, `6`, `7`, `8`), ~replace_na(., 0)) |># I need to replace NA with 0, else I can't sum! According to Olavi who sent the data, NA means 0 and nothing else. Rectangle*year combinations that do not have information about biomass are simply not included in this datamutate(ices_rect =as.factor(ices_rect),Species ="Sprat",biomass_spr =`1`+`2`+`3`+`4`+`5`+`6`+`7`+`8`, IDr =paste(ices_rect, Year, sep ="."),lon =ices.rect(ices_rect)$lon,lat =ices.rect(ices_rect)$lat) |># Make new IDfilter(Year >1992& Year <2020)
mutate: new variable 'sub_div' (character) with 12 unique values and 0% NA
her <- her |>add_utm_columns(ll_names =c("lon", "lat"),utm_crs =32633)# Plot distribution over time in the whole areaplot_map_fc +geom_point(data = spr, aes(X*1000, Y*1000, color = biomass_spr), size =1.5) +facet_wrap(~Year) +scale_color_viridis(trans ="sqrt", name ="Sprat biomass (tonnes)")
# It's because rectangles somehow being in different sub divisions.# I need to group by IDr and summarizespr_sum <- spr |>group_by(IDr) |>summarise(biomass_spr =sum(biomass_spr)) |># Sum abundance within IDrdistinct(IDr, .keep_all =TRUE) |># Remove duplicate IDrmutate(ID_temp = IDr) |># Create temporary IDr that we can use to split in order# to get Year and StatRect back into the summarized dataseparate(ID_temp, c("StatRec", "Year"), sep =4)
group_by: one grouping variable (IDr)
summarise: now 1,418 rows and 2 columns, ungrouped
distinct: no rows removed
mutate: new variable 'ID_temp' (character) with 1,418 unique values and 0% NA
# How many rows per rectangle?spr_sum |>group_by(IDr) |>mutate(n =n()) |>ungroup() |>distinct(n)
group_by: one grouping variable (IDr)
mutate (grouped): new variable 'n' (integer) with one unique value and 0% NA
ungroup: no grouping variables
distinct: removed 1,417 rows (>99%), one row remaining
# A tibble: 1 × 1
n
<int>
1 1
# Now do the same for herringher_sum <- her |>group_by(IDr) |>summarise(biomass_her =sum(biomass_her)) |># Sum abundance within IDrdistinct(IDr, .keep_all =TRUE) |># Remove duplicate IDrmutate(ID_temp = IDr) |># Create temporary IDr that we can use to split in order# to get Year and StatRect back into the summarized dataseparate(ID_temp, c("StatRec", "Year"), sep =4)
group_by: one grouping variable (IDr)
summarise: now 1,418 rows and 2 columns, ungrouped
distinct: no rows removed
mutate: new variable 'ID_temp' (character) with 1,418 unique values and 0% NA
mutate: new variable 'IDr' (character) with 1,769 unique values and 0% NA
Joining with `by = join_by(IDr)`
left_join: added 4 columns (biomass_her, StatRec, Year, biomass_spr)
> rows only in x 131,190
> rows only in y ( 54)
> matched rows 812,238
> =========
> rows total 943,428
# If NA on the rectangle level, take the average values across ices rectangles *within* the sub-division# If no rectangles within the SD, take the annual meancod_pred <- cod_pred |>group_by(year, sub_div) |>mutate(biomass_her_year_sd_mean =mean(biomass_her, na.rm =TRUE),biomass_spr_year_sd_mean =mean(biomass_spr, na.rm =TRUE)) |>ungroup() |>group_by(year) |>mutate(biomass_her_year_mean =mean(biomass_her, na.rm =TRUE),biomass_spr_year_mean =mean(biomass_spr, na.rm =TRUE)) |>ungroup() |>mutate(biomass_her =ifelse(is.na(biomass_her), biomass_her_year_sd_mean, biomass_her),biomass_spr =ifelse(is.na(biomass_spr), biomass_spr_year_sd_mean, biomass_spr)) |>mutate(biomass_her =ifelse(is.na(biomass_her), biomass_her_year_mean, biomass_her),biomass_spr =ifelse(is.na(biomass_spr), biomass_spr_year_mean, biomass_spr)) |> dplyr::select(-biomass_her_year_mean, -biomass_her_year_mean, -biomass_spr_year_sd_mean, -biomass_spr_year_mean)
group_by: 2 grouping variables (year, sub_div)
mutate (grouped): new variable 'biomass_her_year_sd_mean' (double) with 133 unique values and 9% NA
new variable 'biomass_spr_year_sd_mean' (double) with 133 unique values and 9% NA
ungroup: no grouping variables
group_by: one grouping variable (year)
mutate (grouped): new variable 'biomass_her_year_mean' (double) with 28 unique values and 7% NA
new variable 'biomass_spr_year_mean' (double) with 28 unique values and 7% NA
ungroup: no grouping variables
mutate: changed 50,514 values (5%) of 'biomass_her' (50514 fewer NA)
changed 50,514 values (5%) of 'biomass_spr' (50514 fewer NA)
mutate: changed 15,612 values (2%) of 'biomass_her' (15612 fewer NA)
changed 15,612 values (2%) of 'biomass_spr' (15612 fewer NA)
# Test plotplot_map_fc +geom_raster(data = cod_pred |>filter(year <2020), aes(X*1000, Y*1000, fill = biomass_spr)) +facet_wrap(~year) +scale_fill_viridis(trans ="sqrt", name ="Sprat biomass (tonnes)")
# Now add in the sub division values (in case we want to look at that also)biomass_spr_sd <-read_xlsx(paste0(home, "/data/BIAS/N and B per SD 1991-2020.xlsx"),sheet =4) |>mutate(sub_div =ifelse(Sub_Div =="28_2", "28", Sub_Div)) |>filter(sub_div %in%c("24", "25", "26", "27", "28")) |>rename("Year"="ANNUS") |>mutate_at(vars(`AGE1`, `AGE2`, `AGE3`, `AGE4`, `AGE5`, `AGE6`, `AGE7`, `AGE8+`), ~replace_na(., 0)) |># I need to replace NA with 0, else I can't sum! According to Olavi who sent the data, NA means 0 and nothing else. Rectangle*year combinations that do not have information about biomass are simply not included in this datamutate(sub_div =as.factor(sub_div),Species ="Sprat",biomass_spr_sd =`AGE1`+`AGE2`+`AGE3`+`AGE4`+`AGE5`+`AGE6`+`AGE7`+`AGE8+`, # omitting `0+` hereID_sd_year =paste(sub_div, Year, sep =".")) |># Make new ID dplyr::select(biomass_spr_sd, ID_sd_year)
mutate: new variable 'sub_div' (character) with 12 unique values and 0% NA
filter: removed 152 rows (51%), 147 rows remaining
rename: renamed one variable (Year)
mutate_at: changed one value (1%) of 'AGE6' (1 fewer NA)
changed 10 values (7%) of 'AGE7' (10 fewer NA)
changed 13 values (9%) of 'AGE8+' (13 fewer NA)
mutate: converted 'sub_div' from character to factor (0 new NA)
new variable 'Species' (character) with one unique value and 0% NA
new variable 'biomass_spr_sd' (double) with 147 unique values and 0% NA
new variable 'ID_sd_year' (character) with 147 unique values and 0% NA
mutate: new variable 'sub_div' (character) with 12 unique values and 0% NA
filter: removed 152 rows (51%), 147 rows remaining
rename: renamed one variable (Year)
mutate_at: changed 2 values (1%) of 'AGE7' (2 fewer NA)
changed 7 values (5%) of 'AGE8+' (7 fewer NA)
mutate: converted 'sub_div' from character to factor (0 new NA)
new variable 'Species' (character) with one unique value and 0% NA
new variable 'biomass_her_sd' (double) with 147 unique values and 0% NA
new variable 'ID_sd_year' (character) with 147 unique values and 0% NA
# Add in the same id to the pred_gridcod_pred <- cod_pred |>mutate(ID_sd_year =paste(sub_div, year, sep ="."))
mutate: new variable 'ID_sd_year' (character) with 145 unique values and 0% NA
cod_pred <-left_join(cod_pred, biomass_spr_sd)
Joining with `by = join_by(ID_sd_year)`
left_join: added one column (biomass_spr_sd)
> rows only in x 48,144
> rows only in y ( 10)
> matched rows 895,284
> =========
> rows total 943,428
cod_pred <-left_join(cod_pred, biomass_her_sd)
Joining with `by = join_by(ID_sd_year)`
left_join: added one column (biomass_her_sd)
> rows only in x 48,144
> rows only in y ( 10)
> matched rows 895,284
> =========
> rows total 943,428
# Plot. Here I also have NAs. If NA, take the annual averageplot_map_fc +geom_raster(data = cod_pred |>filter(year <2020), aes(X*1000, Y*1000, fill = biomass_spr_sd)) +facet_wrap(~year) +scale_fill_viridis(trans ="sqrt", name ="Sprat biomass (tonnes)")
group_by: one grouping variable (year)
mutate (grouped): new variable 'biomass_her_mean_sd' (double) with 29 unique values and 3% NA
new variable 'biomass_spr_mean_sd' (double) with 29 unique values and 3% NA
ungroup: no grouping variables
mutate: changed 15,612 values (2%) of 'biomass_spr_sd' (15612 fewer NA)
changed 15,612 values (2%) of 'biomass_her_sd' (15612 fewer NA)
# Save plots with "imputed" biomasses# subdivisionplot_map_fc +geom_raster(data = cod_pred |>filter(year <2020), aes(X*1000, Y*1000, fill = biomass_spr_sd)) +facet_wrap(~year) +scale_fill_viridis(trans ="sqrt", name ="Sprat biomass (tonnes)")
plot_map_fc +geom_raster(data = cod_pel |>filter(year <2020),aes(X*1000, Y*1000, fill = cod_her_ovr_sd)) +facet_wrap(~year) +scale_fill_viridis(trans ="sqrt", name ="Cod-Herring overlap (ICES Rect.)",option ="mako",# trim extreme high values to make spatial variation more visiblena.value ="black", limits =c(0, quantile(filter(cod_pel, year <2020)$cod_her_ovr_sd, 0.999)))
